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Abstract 

We investigate the optical response of a system consisting of periodic silver V-grooves interacting 
with quantum emitters. Two surface plasmon-polariton resonances are identified in the reflection 
spectrum of bare silver grooves, with the intensity of one resonance being localized near the bottom 
of the groove and that of the other resonance being distributed throughout the entire groove. 
The linear response of the hybrid silver-emitter system is thoroughly analyzed by considering the 
coupling between surface plasmon polaritons and emitters as the geometry of the grooves and 
the spatial distribution of emitters within the grooves are varied. The nonlinear response of the 
system is also considered by pumping the emitters with a short, high-intensity pulse. By changing 
the duration or the intensity of the pump, the population of emitters in the ground state at the 
end of the pump is varied, and it is found (upon probing with a short pulse) that an increase in 
the fraction of emitters in the ground state corresponds to an increase in Rabi splitting. Spatial 
variations in the ground state population throughout the emitter region are shown to be a result 
of field retardation. 
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I. INTRODUCTION 


Surface plasmon-polaritons (SPPs) are electromagnetic excitations resulting from the cou¬ 
pling of the incident radiation with collective oscillations of conductive electrons near the 
interface between metal and dielectric. Due to the specific dispersion the effective wavelength 
of SPPs is usually significantly shorter than that of the incident field making SPPs highly at¬ 
tractive for various applications 1,2 . These range from electromagnetic energy transport at the 
nanoscale 3 , nano-focusing 4 , utilization of plasmons to achieve lasing beyond the diffraction 
limit 5 , and many others 6 . There is also a growing interest in optics of hybrid nano-materials, 
structures composed of plasmon-polariton sustaining systems (such as nanoparticles, one¬ 
dimensional and two-dimensional periodic arrays, for instance) and molecular aggregates 7 . 
The original interest was geared towards electromagnetic energy exchange between semi- 
classical SPPs and molecular excitons 8 occurring on a femtosecond timescale 9 . Much work 
has been done to investigate linear optical properties of both propagating 10-14 and localized 
SPPs 15-17 coupled to quantum emitters in the strong coupling regime (for a comprehensive 
review of the current state of the field see [7]). Going beyond the linear regime it was re¬ 
cently shown that such systems present an interesting opportunity to actively control light 
at the nanoscale and manipulate optical properties of nano-materials 18,19 . 

Our major point of interest in this work is to scrutinize electromagnetic properties of 
molecular aggregates composed of simple two-level emitters strongly coupled to SPP waves 
supported by plasmonic waveguides. We chose a periodic one-dimensional array of V-grooves 
as an example of such a waveguide. V-grooves are host to various optical phenomena, many 
of which involve the intense, highly localized fields that result from SPP waves. Along these 
lines, the manipulation of channel plasmon-polaritons (CPPs) propagating along metallic 
V-groove channels has been extensively studied 20-22 . The propagation length of CPPs was 
investigated in [20] by filling the grooves with fluorescent microscopic beads and measuring 
the propagation length by microscope, and selection of CPP modes can be accomplished by 
adjusting the shape and size of the groove 21 . Additionally nanofocusing of light, which is of 
extreme practical interest, was demonstrated using tapered V-shaped waveguides 22 . 

Theoretical studies have been conducted using the Green’s Function Integral Equation 
Method 23 , yielding well-defined reflection features that are thought to result from surface 
plasmon polaritons, geometrical resonances, and a Wood’s anomaly 24 . The origins of these 
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features are identified by their response to changes in the geometry of the grooves: period, 
depth, groove angle, and angle of the incident fields. One point of interest is that one 
type of resonance displays intensity that is distributed over the entire groove whereas the 
intensity of another type is localized near the bottom of the groove. In this manuscript we 
numerically investigate how different SPPs interact with molecular aggregates in both linear 
and nonlinear regimes. The latter is considered using pump-probe simulations 25 . 


II. MODEL 


Classical Maxwell’s equations are used to describe the propagation of electromagnetic 
(EM) waves 


dH „ - 

^°~dt = _V X E ’ 

dE - 

e 0~K7 = V X H — J, 
at 


(la) 

(lb) 


where /i 0 and e 0 are the permeability and the permittivity of free space, respectively, E is 
the electric field, H is the magnetic field, and J is the current density either in the metal or 
that due to induced polarization of the emitters as discussed below. 

The finite-difference time-domain (FDTD) method is used to propagate Eqs. (la) and 
(lb) in space and time 26 . All of the results presented herein use a two-dimensional grid in 
which the fields E x , E y and H z are evaluated in the x-y plane and the structure is taken to 
be infinitely long in the .'.--direction, as shown in Fig. 1. 

The system under consideration is open in the y direction and periodic in x. We add 
absorbing boundaries using convolutional perfectly matched layers (CPML) on the top and 
the bottom of the grid as shown in Fig. 1. The left and right sides of the grid are terminated 
with periodic boundary conditions (PBCs) and the incident wave is introduced via a total 
field / scattered field (TFSF) approach 26 . A spatial step of 1.0 nm was selected as results 
are converged for this spacing, which was demonstrated by obtaining the same data using 
a spatial size of 0.5 nm. A time step of dx/(2c) was chosen such that the Courant stability 
condition is satisfied. 

The linear Drude model is used to describe the dispersion of the metal 27 
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FIG. 1 . Schematic diagram of the system. The absorbing boundary conditions are implemented 
as Convolutional Perfectly Matched Layers (CPML), the periodic boundary conditions are applied 
along y-axis, and the incident wave is introduced via a total-field / scattered-field (TFSF) approach. 

where uj p is the plasma frequency, 7 is the phenomenological damping, and e r is the high- 
frequency limit of the dielectric function. For silver, we use the following parameters 28 : 
ujp = 11.59 eV, 7 = 0.2027 eV, and e r = 8.26. In metal, the dynamics of the current density 
J satisfy the following equation 28 : 

<9 7 

— = -7 J + eot 4E. (3) 

In the linear regime, when the frequency response of a system to external EM excitation 
is independent from the incident intensity, one can use a short pulse method to obtain the 
spectrum of the system within a single FDTD run 29 . Under the assumption that only the 
elastic scattering contributes to the spectrum, the reflection is calculated by launching a 
short pulse (of duration r = 0.15 fs and whose form is E 0 sin 2 (y) cos (cot)) and spatially 
integrating the .^-component of the Poynting vector (formed by the cross product of the 
Fourier transformed E- and H-fields; specifically, E X H Z ) along a line of constant y -value on 
the input side. This calculation is performed in the scattered-field region, thereby measuring 
only the reflected fields. 
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The time dynamics of the interaction between the EM field and the molecular aggregate 
is described by the Liouville-von Neumann equation 

= [H,p\-ihfp, (4) 

where p is the single-emitter density matrix, T describes relaxation processes, and H is the 
Hamiltonian describing an emitter with a dipole moment operator jl interacting with an 
electric field 

H = Hq — jl ■ E(t). (5) 

The expectation value of the dipole moment operator is obtained by evaluating Tv(pjl). 
The updated expectation value of the dipole moment is used to calculate the macroscopic 
polarization, P, and then the polarization current, J p , which is subsequently inserted into 
Eq. (lb) 


P = n a (p ), (6a) 

dP 

Jp = i (6b) 

where n a is the volume density of emitters. 

In order to account for all possible local field polarizations we consider quantum emitters 
with three energy levels: an s-type ground state and two degenerate excited p-type states. 
In the basis of angular momentum, wave functions are chosen 29 : |1) = |s), |2) = (| p x ) + 

i \p v ))/V2, 13 ) = (| p x ) ~ i \p y ))/V2. 

The following set of parameters describing a quantum emitter is used in this paper: the 
transition dipole moment is 10 Debye and the radiationless lifetime of the excited state is 1 
ps. The number density and the pure dephasing time are varied. 


III. RESULTS AND DISCUSSION 

The intent of this paper is to scrutinize the optical properties of a periodic system com¬ 
prised of a periodic array of V-grooves in an optically thick silver film (the thickness of 
the film in all simulations is 800 nm) that is optically coupled to quantum emitters. We 
first consider a periodic array of V-grooves without emitters, followed by a hybrid system 
consisting of emitters (all starting in the ground state) added into the grooves, and finally 
this same system with an optical femtosecond pump applied. 
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FIG. 2. Linear optical response of bare V-grooves at normal incidence. Panel (a) shows reflection 
spectrum of bare grooves with 400 nm period, 20 degree groove angle and 200 nm groove depth. 
Panel (b) shows time-averaged intensity in the bare groove when excited by CW plane wave at 1.2 
eV. The intensity is distributed throughout the groove. Panel (c) shows time-averaged intensity 
in the bare groove when excited with CW plane wave at 2.255 eV. The intensity is localized near 
the bottom of the groove. The contour data in panels (b) and (c) is logarithmic and normalized 
to the incident intensity. 


The reflection spectrum of a bare silver grating with a 400 nm period obtained at normal 
incidence is shown in Fig. 2a. Three well-resolved resonances are observed with the energy 
of each depending on the geometrical parameters of the grating. Both of the resonances 
at lower energy are thought to be of a plasmonic, rather than geometrical, nature as each 
disappears when the conductivity of the metal is made infinite (i.e. perfectly reflecting). The 
resonance with the highest energy is a Wood’s anomaly, which is indicated by its wavelength 
corresponding to the groove period. Each of these resonances is observed in [24], where their 
behavior was analyzed by adjusting the geometry of the V-grooves as well as the angle of 
incidence. The incident fields in our simulations are p-polarized with normal incidence, as s- 
polarized waves cannot excite SPP waves. In [24], the field enhancements within the groove 
are substantially smaller for s-polarized waves than p-polarized waves. 

As noted in [24], the intensity of the lower energy resonance is distributed throughout 
the groove whereas that of the higher energy resonance is localized near the bottom (see 
Fig. 2). 

On account of each resonance being well-defined for this geometry, the remainder of this 
work (unless otherwise specified) will focus on V-grooves with a period of 400 nm, a groove 
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angle of 20 degrees, a groove depth of 200 nm, and normally incident fields. 

We now consider the optical response of the system when emitters are added inside the 
grooves. When the emitters are resonant to the structure, normal mode splitting (Rabi 
splitting) is clearly observed (Fig. 3a). The Rabi splitting reaches 325 meV. This amount of 
splitting is considered large for hybrid nanostructures 30 . Either at high emitter densities or a 
large transition dipole moment a third feature appears at or near the emitter resonance as can 
be seen in Fig. 3a near 1.1 eV. This feature, which is not predicted by the coupled-oscillator 
model, has been observed previously in simulations 13,17,25 as well as in experiments 8,31,32 . 
Several observations discussed in [13] suggest that this peak has its origins in SPP enhanced 
emitter-emitter interactions. 

To better understand the physics of this resonance we performed simulations gradually 
varying either groove angle (Fig. 3b) or groove depth (Fig. 3c). Note that the bare SPP 
lines plotted in Fig. 3b and 3c are clearly not linear with groove angle or groove depth 
(see below). This allows us to sweep the SPP resonance through the emitter’s mode. In 
the reflection spectrum we record energy positions of lower and upper branch of the hybrid 
mode. This is carried out at the emitter density of 3 x 10 26 emitters/m 3 (at which the 
third feature is prominent) and the results are shown in Fig. 3b and 3c. The results clearly 
indicate avoided crossing - a unique signature of the strong coupling due to efficient energy 
exchange between the corresponding SPP mode and molecular excitons (in our case these 
are simply two-level emitters). Next, the energy of the intermediate peak (green triangles 
in Fig. 3b and c) does not deviate appreciably from the emitter resonance even as the 
SPP resonance is tuned. Furthermore, this peak merges with the upper polariton as the 
thickness of a spacer layer between the grating and emitters increases. The dipole coupling 
between the emitters themselves is therefore suspected given the fall-off of this peak as the 
SPP field at the emitters’ location decreases. Finally, a simulation was run in which the 
entire region containing emitters was replaced by a single two-level system (essentially a 
spatially distributed single emitter with a very large dipole moment), thereby eliminating 
any possible interaction between the emitters. In this case, the third peak disappears even 
under extremely high coupling conditions. This confirms earlier findings 13,17 which suggested 
that the intermediate resonance located in the middle of the Rabi splitting corresponds to 
dipole-dipole interactions between emitters greatly enhanced by the SPP mode. 

The dashed lines in Fig. 3b and 3c are calculations of the values of the upper and lower 
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polaritons using the coupled oscillator model. The Hamiltonian of the coupled system is 


E m A 


( 7 ) 


A E p i 

E p i is the SPP energy of the bare metallic grooves (obtained from simulations), E m is the 
transition energy of uncoupled emitters, and 2A is the minimum Rabi splitting value. The 
eigenvalues obtained are 

(Epi + E m ) ± yf (Epi + T'm ) 2 — 4A 2 


Eu,l — 


( 8 ) 


Eu,l are the energies of the upper or lower polaritons. In Fig. 3b, the comparison to the 
coupled oscillator model is close whereas in Fig. 3c it deviates somewhat at larger groove 
depths. In particular, we note that both the upper and lower polaritons appear to be 
’’pinched” toward each other in Fig. 3c. We offer two possible explanations for this. First, a 
large groove depth requires the incident fields to traverse a larger length of emitters. More 
absorption takes place than for a shallower grove, leading to decreased excitation of SPP 
waves near the bottom and therefore less coupling. Second, the third feature may interact 
with the upper and lower polaritons at greater groove depths in such a way as to reduce the 
coupling. 

While it is clear that many of the emitters are coupled to the SPPs, some may remain 
coupled only to the incident field 33 and thus act as an absorbing layer. To better understand 
the overall optical coupling in spatially distributed inhomogeneous hybrid systems, two 
different arrangements of emitters are simulated: ’’full” grooves and ’’hollowed” grooves, as 
shown in Fig. 4a. In the ’’full” grooves, the region of the grooves occupied by emitters is 
completely full of emitters up to a given height, whereas in ’’hollowed” grooves, the region 
of emitters extends out sideways from either side of the groove by a fixed width (referred to 
as ’’width of emitters” from here on). The grooves are illuminated with CW fields at 1.2 eV 
(uniform intensity, Fig. 2b) or 2.255 eV (localized intensity, Fig. 2c). 

The Rabi splitting is observed as the area occupied by the emitters is varied. Fig. 4b 
shows that, for the distributed resonance, there is relatively little difference in Rabi splitting 
when the same amount of area is occupied by the emitters for either full or hollowed grooves, 
which demonstrates that the entire volume of emitters is indeed coupled to the SPP waves. 
Given that coupling strength depends on field strength, this is reasonable in light of the 








Groove Angle (Degrees) Groove Depth (nm) 

FIG. 3. Optics of periodic V-grooves coupled to quantum emitters. Panel (a) shows the reflection 
for bare grooves (dashed line) and the reflection for grooves with emitters is shown as a solid 
line (with emitters resonant at 1.2 eV). Panel (b) shows the upper polariton (red squares), lower 
polariton (blue circles), and third resonant mode (green triangles) as a function of the groove angle 
to sweep the SPP resonance through the emitter resonance. Panel (c) shows the same modes as 
in panel (b) but as functions of the groove depth to sweep the SPP resonance through the emitter 
resonance. The horizontal black lines in panels (b) and (c) represent the fixed emitter resonance, 
the curved black lines represent the SPP energy of the bare grooves for the given geometry, and 
the dashed lines are values predicted by the coupled oscillator model. The density of emitters is 
3 x 10 26 emitters/m 3 . 


uniform distribution of intensity in the bare groove at this frequency; the same Rabi splitting 
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Emitter Area (nm 2 ) Emitter Area (nm 2 ) 

FIG. 4. Spatially dependent coupling. The density of emitters is 10 26 emitters/m 3 , the transition 
energy is 1.2 eV, and the pure dephasing time is 400 fs. As the height of the molecular aggregate 
is increased, the coupling for the distributed resonance increases continuously whereas that for the 
localized resonance levels off quickly. Panel (a) shows schematics of’’full” grooves versus ’’hollowed” 
grooves. Panel (b) shows the Rabi splitting as a function of the area occupied by the emitters for 
the more spatially distributed 1.2 eV resonance, where blue circles indicate full grooves, red squares 
indicate hollowed grooves (width of emitter region is lOnm) and green triangles indicate hollowed 
grooves (width of emitter region is 15 nm). Panel (c) is the same as panel (b) except the data is 
shown for the more localized 2.255 eV resonance. The pure dephasing time is 600 fs. 


is achieved regardless of where a given area of emitters is placed in the groove. Fig. 4c shows 
that, for the localized resonance, a large increase in Rabi splitting occurs as the area is 
increased from the bottom, but then it levels off to a constant value well before the emitters 
reach the top of the groove. This is exactly what is expected as the strong fields near the 
bottom of the groove give the strongest coupling. The constant value of Rabi splitting is 
smaller for the lowest width of emitters (10 nm) and this occurs because the fields above 
the bottom of the groove are weaker but not zero (hence a larger hollow region loses some of 
the emitters coupled to those weaker fields). A small jump occurs at the end of each graph 
in Fig. 4c and this is due to the increased fields at the sharp corners of the groove. This 
jump occurs at smaller areas for more hollowed grooves as they occupy less area when they 
extend to the top of the groove. 

All simulations discussed above begin with all of the emitters in the ground state. One can 
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FIG. 5. Pump-probe dynamics. Panel (a) shows Rabi splitting as a function of pump amplitude. 
Each value is obtained by launching a short, low-intensity probe pulse immediately after the pump. 
Panel (b) shows ground state population, p \\, as a function of pump amplitude. This value is 
averaged over the entire region occupied by emitters and obtained at the end of the pump. The 
density of emitters is 10 26 emitters/m 3 , the pure dephasing time is 400 fs, the transition energy is 
1.2 eV, and the duration of the pump is 30 fs. 

pump the system by sending in a high-intensity pulse, thereby inducing Rabi oscillations 
in the emitters. Below the time dynamics of a pumped hybrid system is discussed. In 
particular, we observe how Rabi splitting depends on the ground state population at the 
end of the pump. 

First, the system is pumped with a 30 fs pulse, and it is subsequently probed with a 
short, low-intensity pulse as it was in the linear regime. The Rabi splitting is obtained from 
the reflection spectrum by calculating the difference in resonant energies for the upper and 
lower polaritons. Fig. 5a shows Rabi splitting as a function of pump amplitude at the end 
of the pump. Fig. 5b shows the ground state population, averaged over the entire region 
of emitters at the end of the pump, as a function of pump amplitude. Notice how the plot 
in Fig. 5b oscillates; this is because the area under the pump pulse (which depends on 
both pump amplitude and duration) determines the number of Rabi cycles 34 . Thus different 
pump amplitudes generate greater or fewer Rabi oscillations, leading to different values of 
the ground state population at the end of the pump. 

Fig. 5 indicates that the amount of Rabi splitting depends on the ground state population 
of the emitters, pn. We see that a large excited population of emitters in the entire groove 
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will reduce or eliminate the Rabi splitting whereas a smaller excited population will yield 
larger Rabi splitting. We speculate that the coupling may take on a different character in 
the non-linear regime, perhaps to the extent that Rabi splitting is not the only indicator of 
coupling strength. Further investigation of the optical properties of excited emitters coupled 
to plasmons is clearly warranted. 

Although the changes of pn are in step with changes in Rabi splitting, large decreases 
of pn do not always give a correspondingly large decrease in Rabi splitting and this can be 
understood in terms of the non-uniformity of pn throughout the groove. In particular, the 
pumping fields arrive at deeper parts of the groove later than they arrive at the top and the 
Rabi oscillations are not perfectly in phase with one another along the length of the groove 
due to retardation. Additionally, the SPP fields are inhomogeneous throughout the groove. 
Thus, the spatially averaged value of pn after the pump gives a good, though not ideal, 
indication of the subsequent coupling to SPP fields. 

Fig. 6a shows a spatial distribution of pn at the end of 70 fs long pump. A clear 
strong spatial variation of the ground state population is seen. In fact, the variations of the 
ground state population are oscillations whose wavelength varies somewhat over the region 
of emitters. One might surmise that the wavelength of these oscillations is on the order of 
that of the pump, but it is actually significantly smaller. To investigate this further, we 
numerically solve the Schrodinger equation for a one-dimensional region of two-level atoms 
(finite along, say, the 2 -axis while infinite along two others) subject to excitation from a 
pump. 

Because we consider a one-dimensional region, retardation effects must be included by 
using the retarded time t — z/v, where v is the speed of light in the medium in the expression 
for the pulse. The retardation effects are revealed to be the cause of the spatial modulations 
of the ground state population, as shown in Fig. 6b: the temporal oscillations of the ground 
state probability are similar between adjacent spatial points, but slightly shifted. Hence at 
a given time, adjacent points have slightly different values of ground state probability. To 
further elucidate this idea, we generated larger phase shifts in the temporal oscillations be¬ 
tween nearby points by adjusting two parameters: the pump amplitude and the propagation 
velocity of light in the emitter region. 

For a larger pump amplitude, each emitter will undergo Rabi flopping more rapidly in 
time. Thus a given phase shift between adjacent locations will lead to a larger shift in the 
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FIG. 6. Spatial modulations of the molecular ground state. Each run consists of a 70 fs pump 
applied to the hybrid V-grooves system. The simulations shown in panels (a), (c), and (d) are 
run with FDTD and the Liouville-von Neumann equation whereas that in panel (b) is run by 
numerically integrating the Schrodinger equation for a two-level atom. For panels (a), (c), and (d) 
the density is 10 26 emitters/m 3 , the pump amplitude is 7 x 10 8 V/m, and the dephasing time is 400 
fs. Panel (a) shows the fraction of emitters in the ground state. The characteristic length of spatial 
modulations is much less than the pump wavelength of 1033 nm. Panel (b) shows oscillations of 
the ground state population for two emitters separated by 50 nm. Panel (c) shows that the group 
velocity at the transition frequency is less than c. Panel (d) shows the ground state population as 
a function of coordinate for the system pumped on resonance (1.2 eV, solid line) and slightly off 
resonance (1.15 eV, dashed line). 


ground state population between those two locations. This was observed in our simulations: 
in general, increasing the pump amplitude leads to an increase in the number of spatial 
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oscillations of the ground state population over the region of emitters. For a slower group 
velocity, the pump takes longer to reach an adjacent point, causing the temporal oscillations 
between two adjacent points to acquire a larger phase difference. Our simulations allow 
for the adjustment of the group velocity of light in the emitter region, and we see that a 
decrease in the group velocity results in a greater number of spatial oscillations relative to 
a larger speed of light. 


Fig. 6c shows that group velocity is decreased (in fact, negative since most emitters are 
inverted) at the transition frequency of 1.2 eV, and it is higher away from resonance. A neg¬ 
ative group velocity is possible in inverted systems and has been experimentally observed 35 . 
It is stated that this occurs when different frequency components of a pulse interfere with 
one another (in a region of anomalous dispersion) in such a way as to cause a resonant pulse 
to be advanced relative to a non-resonant pulse traveling at c. We are assured by [36] that 
this does not violate causality and that it occurs because the early parts of the pulse are 
reshaped to resemble the later pulse. Furthermore, [37] points out that in passing through a 
medium with negative group velocity, the information transmitted by the pulse front suffers 
a positive and causal delay. Our system was pumped (in separate runs) both on and off 
of resonance and Fig. 6d shows that the wavelength of the spatial oscillations for the off 
resonance pump is indeed increased relative to the on resonance pump on account of the 
latter having a slower group velocity. 


All of these results indicate that retardation effects give rise to spatial oscillations in the 
ground state population of the emitters contained within the groove. The larger fields of 
the surface plasmons and the slower group velocity within the emitter region surely enhance 
these spatial oscillations. It may be possible to use this phenomenon to perform a new 
type of optical design using femtosecond pulses as a tool to craft hybrid systems. Highly 
inhomogeneous spatial modulations of molecules lead to the modified refractive index, which 
is appreciably anisotropic. One can envision an exciting opportunity for a new research 
direction, in which both the geometry of metal nanostructures and pump pulses govern the 
refractive index of the system. Furthermore one can apply optimization techniques such as 
genetic algorithms 38 , for instance, to design materials with desired optical properties. 
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IV. CONCLUSION 


The optical properties of the bare and hybrid V-groove systems have been explored un¬ 
der several different circumstances. Two SPP resonances with entirely different spatial 
distributions were shown. Optical coupling between quantum emitters and SPP waves was 
thoroughly characterized by simulating different groove geometries and spatial configura¬ 
tions of emitters’ distributions within the grooves. The value of Rabi splitting was shown 
to vary with the pumping intensity and in accord with the ground state population. An 
explanation of this awaits a more developed understanding of coupling in the non-linear 
regime. Spatial oscillations of the ground state population of emitters within the groove 
are shown to be the result of field retardation. This work puts forth several suggestions for 
experimental investigation of coupling in hybrid systems, most notably an investigation of 
the time dynamics of coupling when the system is pumped with different intensities. 
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